Análisis peces kelp forest

suppressPackageStartupMessages({
  library(vegan)
  library(here)
  library(magrittr)
  library(cowplot)
  library(sf)
  library(tmap)
  library(tidyverse)
})

kelp <- read.csv(here("data", "Peces_KelpForest_2011-2013.csv"),
                 stringsAsFactors = F)

kelp11 <- kelp %>% 
  filter(year == 2011)

Pregunta 1: Existen diferencias en la composición de la comunidad?

Para datos 2011

Sitios con transectos Fondo y Media Agua

kelp11 %>%
  group_by(location, site, level, transect) %>%
  tally() %>% 
  ungroup() %>%
  select(-n) %>% 
  group_by(location, site, level) %>% 
  tally() %>% 
  knitr::kable(caption = "Numero de transectos por localidad, sitio (N-S) y nivel (F-MA), incluyendo registros de 'fuera de transecto'.")

Numero de transectos por localidad, sitio (N-S) y nivel (F-MA), incluyendo registros de ‘fuera de transecto’.

location	site	level	n
ASA	N	Bottom	4
ASA	N	Midwater	4
ASA	S	Bottom	4
ASA	S	Midwater	3
BMA	N	Bottom	4
BMA	N	Midwater	4
BMA	S	Bottom	4
BMA	S	Midwater	4
CKE	N	Bottom	9
CLO	N	Bottom	3
CLO	N	Midwater	2
CLO	S	Bottom	4
CLO	S	Midwater	4
COL	N	Bottom	1
COL	S	Bottom	2
ERE	N	Bottom	4
ERE	N	Midwater	4
ERE	S	Bottom	4
ERE	S	Midwater	2
ERO	N	Bottom	4
ERO	N	Midwater	4
ERO	S	Bottom	4
ERO	S	Midwater	4
ISME	N	Bottom	4
ISME	N	Midwater	4
ISME	S	Bottom	4
ISME	S	Midwater	4
ITSP	N	Bottom	4
ITSP	N	Midwater	4
ITSP	S	Bottom	4
ITSP	S	Midwater	4
PCH	S	Bottom	4
PCH	S	Midwater	4
RET	N	Bottom	4
RET	N	Midwater	4
RET	S	Bottom	4
RET	S	Midwater	3
SMI	N	Bottom	4
SMI	N	Midwater	4
SMI	S	Bottom	4
SMI	S	Midwater	4
SSI	N	Bottom	4
SSI	N	Midwater	4
SSI	S	Bottom	4
SSI	S	Midwater	4
STO	N	Bottom	3
STO	N	Midwater	3
STO	S	Bottom	3
STO	S	Midwater	3
VTR	N	Bottom	3
VTR	N	Midwater	2
VTR	S	Bottom	1
VTR	S	Midwater	1

De la tabla anterior vemos que Campo Kenedy, Colonet, Campo Lopez, Eréndira y Vaye Tranquilo no tienen suficientes muestras (mínimo 3 transectos por nivel por zona), por lo que las excluimos en el siguiente análisis.

Todas las matrices de distancias son calculadas por distancia de Bray-Curtis.

Diferencias presencia / ausencia

Necesitamos una matriz de presencia / ausencia por especie para correr el ANOSIM entre F y MA

localidades_excluidas <- c("CLO", "COL", "ERE", "VTR", "CKE", "PCH")

data_test1 <- kelp11 %>% 
  filter(!location %in% localidades_excluidas) %>% 
  group_by(location, site, level, latitude, longitude, genusspecies) %>% 
  summarize(abundance = sum(abundance, na.rm = T)) %>% 
  ungroup() %>% 
  mutate(abundance = 1) %>% 
  spread(genusspecies, abundance, fill = 0)

data_test1_groups <- data_test1 %>% 
  select(location, site, level) %>% 
  mutate(loc_site = paste(location, site, sep = "-"))

data_test1_samples <- data_test1 %>% 
  select(-c(location, site, level, latitude, longitude)) %>% 
  vegdist(method = "bray")

Habiendo filtrado las localidades que no tenian suficientes muestras, obtenemos la siguiente lista de localidades: ASA, BMA, ERO, ISME, ITSP, RET, SMI, SSI, STO.

Los datos de arriba estan listos para correr los ANOSIMs. El objeto data_test1_groups tiene los factores o tratamientos correspondientes a cada combinacion localidad - sitio - nivel. Los datos tienen la suma a traves de los 4 transectos de cada localidad - sitio - nivel, y son convertidos a presencia / ausencia.

A continuación los resultados de ADONIS (ANOSIMS con ANOVAS de matrices de distancia), que son más robustos para anidaciones (diferente a dos vías).

ANOSIM (ADONIS) entre niveles F y MA

Anidados a nivel localidad-sitio

Reconociendo que las diferencias entre Fondo y MA pueden ser generadas por procesos especificos a cada localidad, corremos el ANOSIM anidando los transectos F y MA dentro de localidad - sitio. En este caso, el ANOSIM anidado indica que dentro de cada localidad - sitio, existen diferencias entre el fondo y media agua.

perm <- how(nperm = 999)
setBlocks(perm) <- data_test1_groups$loc_site

set.seed(43)
adonis(formula = data_test1_samples ~ level, data = data_test1_groups, permutations = perm)

## 
## Call:
## adonis(formula = data_test1_samples ~ level, data = data_test1_groups,      permutations = perm) 
## 
## Blocks:  data_test1_groups$loc_site 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)   
## level      1    0.1391 0.13906  1.2869 0.03647   0.01 **
## Residuals 34    3.6740 0.10806         0.96353          
## Total     35    3.8131                 1.00000          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anidando localidad

La comparacion similar anidando a nivel de localidad (ignorando las diferencias entre N y S) indica que dentro de cada localidad, independiente si es sitio N o S, también hay diferencias entre F y MA.

setBlocks(perm) <- data_test1_groups$location

set.seed(43)
(test1a <- adonis(formula = data_test1_samples ~ level, data = data_test1_groups, permutations = perm))

## 
## Call:
## adonis(formula = data_test1_samples ~ level, data = data_test1_groups,      permutations = perm) 
## 
## Blocks:  data_test1_groups$location 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)  
## level      1    0.1391 0.13906  1.2869 0.03647  0.034 *
## Residuals 34    3.6740 0.10806         0.96353         
## Total     35    3.8131                 1.00000         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Sin anidar

Sin embargo, podemos correr tambien la comparacion de todos los fondos contra todos los MA, ignorando la variabilidad explicada por cada localidad - sitio. En este caso, las comunidades (determinadas por presencia / ausencia) no son suficientemente diferentes.

set.seed(43)
adonis(formula = data_test1_samples ~ level, data = data_test1_groups)

## 
## Call:
## adonis(formula = data_test1_samples ~ level, data = data_test1_groups) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)
## level      1    0.1391 0.13906  1.2869 0.03647  0.268
## Residuals 34    3.6740 0.10806         0.96353       
## Total     35    3.8131                 1.00000

ANOSIM entre Sitio N y S

Ahora debenos comparar entre sitios, anidando a nivel de localidad. Es decir, comparar dentro de cada localidad los sitios N contra S, sin tomar en cuenta la variación que hay entre F y MA. Los resultados indican que para cada localidad, N y S (juntando F y MA) son suficientemente similares como para considerar cada localidad como un grupo homogéneo. No hacemos la comparación N - S sin anidar, pues N - S solamente se definen así dentro de cada localidad, pero no representan fuentes sistemáticas de variación dentro de cada sitio. La historia sería diferente si fuera Expuesto - Protegido, pero eso es en los datos de 2013, no en estos.

setBlocks(perm) <- data_test1_groups$location

set.seed(43)
adonis(formula = data_test1_samples ~ site, data = data_test1_groups, permutations = perm)

## 
## Call:
## adonis(formula = data_test1_samples ~ site, data = data_test1_groups,      permutations = perm) 
## 
## Blocks:  data_test1_groups$location 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model     R2 Pr(>F)
## site       1    0.0968 0.09684 0.88599 0.0254  0.195
## Residuals 34    3.7162 0.10930         0.9746       
## Total     35    3.8131                 1.0000

ANOSIM entre localidades

Finalmente, comparamos entre localidades. En este caso no hay anidacion, pues la localidad es el elemento jerárquico más alto. Los resultados de ADONIS indican diferencias significativas entre sitios.

set.seed(43)
adonis(formula = data_test1_samples ~ location, data = data_test1_groups)

## 
## Call:
## adonis(formula = data_test1_samples ~ location, data = data_test1_groups) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs  MeanSqs F.Model      R2 Pr(>F)    
## location   8    1.9286 0.241076  3.4541 0.50579  0.001 ***
## Residuals 27    1.8845 0.069795         0.49421           
## Total     35    3.8131                  1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Comparación de densidades promediadas

Necesitamos obtener la densidad promedio de cada especie a nivel de localidad - sitio - nivel. En este caso es necesario filtrar los datos de “fuera de transecto”, o transecto 0. Las densidades promedio son transformadas por raiz cuadrada.

data_test2 <- kelp11 %>% 
  filter(transect > 0,
         !location %in% localidades_excluidas) %>% 
  group_by(location, site, level, transect, genusspecies) %>% 
  summarize(abundance = sum(abundance)) %>% 
  ungroup() %>% 
  group_by(location, site, level, genusspecies) %>%
  summarize(abundance = mean(abundance, na.rm = T)) %>%
  ungroup() %>%
  mutate(abundance = sqrt(abundance)) %>% 
  spread(genusspecies, abundance, fill = 0)

data_test2_groups <- data_test2 %>% 
  select(location, site, level) %>% 
  mutate(loc_site = paste(location, site, sep = "-"))

data_test2_samples <- data_test2 %>% 
  select(-c(location, site, level)) %>% 
  vegdist(method = "bray")

nMDS para detectar muestras “problematicas”

Primero con los promedios por localidad - sitio - nivel (i.e. 4 puntos por localidad)

set.seed(43)
mds <- metaMDS(data_test2_samples, trace = F)

stress <- paste("2D Stress =", formatC(mds$grstress, digits = 4, format = "f"))

p1 <- cbind(data_test2_groups, scores(mds)) %>% 
  ggplot(aes(x = NMDS1, y = NMDS2, sitio = site)) +
  geom_point(size = 4, aes(color = location, shape = level)) +
  coord_equal() +
  scale_color_brewer(palette = "Paired") +
  annotate(geom = "text", x = 0, y = 0.5, label = stress)

p1

Version animada:

plotly::ggplotly(p1)

El nMDS muestra claramente diferencias en la densidades entre F y MA.

Podemos corroborar haciendo un nMDS con todos los transectos (sin promediar). En este caso, cada localidad deberia de tener al rededor de 12 puntos (3 NF, 3NM, 3SF, 3SM). Las siguientes figuras es exactamente la misma informacion, pero con diferente representacion. La primera solamente tiene informacion visual de nivel (F vs MA), la segund aincluye tambien color segun localidad.

data_test2 <- kelp11 %>% 
  filter(transect > 0,
         !location %in% localidades_excluidas) %>% 
  group_by(location, site, level, transect, genusspecies) %>% 
  summarize(abundance = sum(abundance)) %>% 
  ungroup() %>% 
  mutate(abundance = sqrt(abundance)) %>% 
  spread(genusspecies, abundance, fill = 0)

data_test2_groups <- data_test2 %>% 
  select(location, site, level) %>% 
  mutate(loc_site = paste(location, site, sep = "-"))

data_test2_samples <- data_test2 %>% 
  select(-c(location, site, level)) %>% 
  vegdist(method = "bray")

set.seed(43)
mds <- metaMDS(data_test2_samples, trace = F)

stress <- paste("2D Stress =", formatC(mds$grstress, digits = 4, format = "f"))

p1 <- cbind(data_test2_groups, scores(mds)) %>% 
  ggplot(aes(x = NMDS1, y = NMDS2, sitio = site)) +
  geom_point(size = 4, aes(shape = level, color = level)) +
  coord_equal() +
  scale_color_brewer(palette = "Paired") +
  annotate(geom = "text", x = 0, y = 0.5, label = stress)

p2 <- cbind(data_test2_groups, scores(mds)) %>% 
  ggplot(aes(x = NMDS1, y = NMDS2, sitio = site)) +
  geom_point(size = 4, aes(color = location, shape = level)) +
  coord_equal() +
  scale_color_brewer(palette = "Paired") +
  annotate(geom = "text", x = 0, y = 0.5, label = stress)

plot_grid(p1, p2, ncol = 1, labels = "AUTO")

ANOSIM 1 via F vs MA

El ANOSIM confirma las diferencia sobservadas en el nMDS de arriba. Este anosim compara las diferencias entre fondo y media agua, anidando por localidad y sitio

perm <- how(nperm = 999)
setBlocks(perm) <- data_test2_groups$loc_site

set.seed(43)
adonis(formula = data_test2_samples ~ level, data = data_test2_groups, permutations = perm)

## 
## Call:
## adonis(formula = data_test2_samples ~ level, data = data_test2_groups,      permutations = perm) 
## 
## Blocks:  data_test2_groups$loc_site 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##            Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## level       1    3.2863  3.2863  17.517 0.14656  0.001 ***
## Residuals 102   19.1356  0.1876         0.85344           
## Total     103   22.4219                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ANOSIM 2 vias

solo F sitio y localidad

Usando solamente los datos de FONDO, corremos un ANOSIM de dos vias comparando entre localidad y sitio. En este caso identificamos que no hay diferencias entre sitio, pero si entre localidades.

# Filtramos los datos para tener solamente FONDO
data_test2_groups_B <- data_test2 %>% 
  select(location, site, level) %>% 
  filter(level == "Bottom") %>% 
  mutate(loc_site = paste(location, site, sep = "-"))

data_test2_samples_B <- data_test2 %>% 
  filter(level == "Bottom") %>%
  select(-c(location, site, level)) %>% 
  vegdist(method = "bray")

set.seed(43)
adonis(formula = data_test2_samples_B ~ site + location, data = data_test2_groups_B)

## 
## Call:
## adonis(formula = data_test2_samples_B ~ site + location, data = data_test2_groups_B) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## site       1    0.2028 0.20279  1.7547 0.02111  0.086 .  
## location   8    4.3178 0.53973  4.6702 0.44951  0.001 ***
## Residuals 44    5.0850 0.11557         0.52938           
## Total     53    9.6056                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

solo M.A. sitio y localidad

Usando solamente los datos de MEDIA AGUA, corremos un ANOSIM de dos vias comparando entre localidad y sitio. En este caso, las diferencias entre localidades son significativas, pero no las diferencias entre sitio.

# Filtramos los datos para tener solamente MEDIA AGUA
data_test2_groups_M <- data_test2 %>% 
  select(location, site, level) %>% 
  filter(level == "Midwater") %>% 
  mutate(loc_site = paste(location, site, sep = "-"))

data_test2_samples_M <- data_test2 %>% 
  filter(level == "Midwater") %>%
  select(-c(location, site, level)) %>% 
  vegdist(method = "bray")

set.seed(43)
adonis(formula = data_test2_samples_M ~ site + location, data = data_test2_groups_M)

## 
## Call:
## adonis(formula = data_test2_samples_M ~ site + location, data = data_test2_groups_M) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## site       1    0.2619 0.26189  1.7431 0.02748  0.100 .  
## location   8    3.2584 0.40731  2.7110 0.34192  0.001 ***
## Residuals 40    6.0096 0.15024         0.63060           
## Total     49    9.5300                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Nivel y Localidad

Ahora comparamos la composicion de la comunidad entre niveles y entre localidades. Esto nos indica que hay diferencias entre ambos.

set.seed(43)
adonis(formula = data_test2_samples ~ level + location, data = data_test2_groups)

## 
## Call:
## adonis(formula = data_test2_samples ~ level + location, data = data_test2_groups) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##            Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## level       1    3.2863  3.2863 21.8719 0.14656  0.001 ***
## location    8    5.0121  0.6265  4.1698 0.22354  0.001 ***
## Residuals  94   14.1235  0.1502         0.62990           
## Total     103   22.4219                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pregunta 2: Es posible identificar una transición biogeografica en los peces de kelp?

Pregunta 3: Que especies explican mejor las diferencias?

Contestamos ambas preguntas con el IVB y SIMPER

IVB

source('~/GitHub/bvi/bvi.R')
source('~/GitHub/bvi/bvi_plot.R')
source('~/GitHub/bvi/bvi_col.R')
source('~/GitHub/bvi/bvi_boxplot.R')

Para el IVB, usamos unicamente los datos del 2011 y con las mismas localidades que hemos usado hasta ahora. (Al incluir todas las localidades del 2011 los resultados son los mismos). Calculamos las densidades promedio por localidad y especie. Cuando una especie no estaba presente, recibe densidad de 0 (Esto se puede cambiar, pero los resultados son los mismos).

site_data <- kelp %>% 
  group_by(year, location, site, latitude, longitude) %>% 
  tally()

bvi_results <- filter(kelp, transect > 0) %>%
  filter(year == 2011,
         !location %in% localidades_excluidas) %>% 
  group_by(location, site, level, transect, genusspecies) %>%
  summarize(n = sum(abundance)) %>%
  spread(location, n, fill = 0) %>%
  gather(location, n,-c(site, level, transect, genusspecies)) %>%
  group_by(location, genusspecies) %>%
  summarize(n = mean(n)) %>%
  ungroup() %>% 
  select(location, genusspecies, n) %>%
  spread(location, n) %>%
  rename(Spp = genusspecies) %>% 
  bvi(sum = F, others = T)

knitr::kable(bvi_results)

Spp	ASA	BMA	ERO	ISME	ITSP	RET	SMI	SSI	STO	BVI	rBVI
Oxyjulis californica	17	16	11	17	17	15	17	16	17	143	10.206995
Semicossyphus pulcher	16	15	12	16	16	14	14	14	15	132	9.421841
Chromis punctipinnis	12	17	16	6	15	0	16	15	10	107	7.637402
Embiotoca jacksoni	0	14	7	0	13	17	15	12	11	89	6.352605
Paralabrax clathratus	15	8	5	15	9	7	7	7	7	80	5.710207
Embiotoca lateralis	0	0	9	14	11	16	0	6	16	72	5.139186
Oxylebius pictus	14	7	0	4	10	13	9	2	7	66	4.710921
Sebastes atrovirens	6	10	10	2	8	9	0	5	12	62	4.425410
Hypsypops rubicundus	0	9	0	13	7	7	8	9	0	53	3.783012
Rhacochilus vacca	0	5	13	9	0	0	0	4	13	44	3.140614
Hypsurus caryi	7	1	6	0	12	0	0	0	9	35	2.498216
Sebastes caurinus	9	4	0	6	0	2	0	0	8	29	2.069950
Brachyistius frenatus	10	0	15	0	0	0	3	0	0	28	1.998572
Rhinogobiops nicholsii	0	1	0	0	14	0	0	10	2	27	1.927195
Sebastes serranoides	0	3	0	0	0	3	11	8	0	25	1.784440
Hyperprosopon anale	0	0	17	0	0	7	0	0	0	24	1.713062
Girella nigricans	0	12	0	10	0	0	0	1	0	23	1.641684
Others	32	32	34	42	27	56	59	44	36	362	25.838687

Visualizar IVB

En BOXPLOT

Cada punto es el “score” de la especie para cada localidad (tabla de arriba).

p1 <- bvi_results %>% 
    select(-c(BVI, rBVI)) %>%
    gather(Sample, Score, -1) %>% 
    set_colnames(value = c("Spp", "Sample", "Score")) %>% 
    transform(Spp = reorder(Spp, Score)) %>% 
    ggplot(aes(x = Spp, y = Score)) +
    geom_boxplot(outlier.shape = NULL,
                 outlier.alpha = 0,
                 fill = "gray",
                 alpha = 0.5) +
    geom_point(aes(fill = Sample), pch = 21, color = "black", size = 2) +
    theme_bw() +
    coord_flip() +
    labs(x = "Spp", y = "Score") +
  scale_fill_brewer(palette = "Paired")

p1

Version animada

plotly::ggplotly(p1)

Podemos visualizar los resultados del IBV (y buscar transiciones en las especies). En este caso, las localidades estan arregladas latitudinalmente, y las especies segun su importancia determinada por el IVB. Las dos figuras son iguales, pero con “smoothing” en la superficie. Sobre todo en el centro, se puede ver como algunas especies desaparecen (al movernos horizontalmente) y otras aparecen. Las de hasta arriba son las constantes, las de hasta abajo son las raras.

bvi_results_plot <- bvi_results %>% 
  filter(!Spp == "Others") %>% 
  gather(Location, Score, -c(BVI, rBVI, Spp)) %>% 
  select(Location, Spp, Score, BVI, rBVI) %>% 
  left_join(site_data, by = c("Location" = "location")) %>% 
  transform(Location = reorder(Location, -latitude)) %>% 
  transform(Spp = reorder(Spp, rBVI))

ggplot(bvi_results_plot, aes(x = Spp, y = Location, fill = Score)) +
  geom_raster() +
  coord_flip() +
  scale_fill_gradientn(colours = colorRamps::matlab.like(20))

ggplot(bvi_results_plot, aes(x = Spp, y = Location, fill = Score)) +
  geom_raster(interpolate = T) +
  coord_flip() +
  scale_fill_gradientn(colours = colorRamps::matlab.like(20))

SIMPER

Ahora hacemos el SIMPER con una aproximacion similar. En este caso, las densidades las transformamos por raiz cuadrada y corremos 999 permutaciones. Al graficar, presentamos solamente las especies que contribuyen al 95% de las disimilitudes. En el eje x presentamos la disimilitud promedio explicada por cada especie.

simper_data <- filter(kelp, transect > 0) %>%
  filter(year == 2011,
         !location %in% localidades_excluidas) %>% 
  group_by(location, site, level, transect, genusspecies) %>%
  summarize(n = sum(abundance)) %>%
  ungroup() %>% 
  mutate(n = sqrt(n)) %>% 
  spread(genusspecies, n, fill = 0)


comm <- simper_data %>% 
  select(-c(location, site, level, transect)) %>% 
  as.matrix()

sim <- simper(comm = comm,
              group = simper_data$location,
              permutations = 999,
              parallel = 3)

cusums <- tibble(species = sim$ASA_BMA$species) %>% 
  cbind(map_df(sim, "cusum")) %>% 
  gather(pair, cumsum, -species)

overall <- map_df(sim, "overall") %>% 
  gather(pair, overall)

p1 <- tibble(species = sim$ASA_BMA$species) %>% 
  cbind(map_df(sim, "average")) %>% 
  gather(pair, average, -species) %>%
  left_join(cusums, by = c("species", "pair")) %>% 
  left_join(overall, by = "pair") %>% 
  filter(cumsum <= 0.95) %>%
  transform(species = reorder(species, average)) %>% 
  mutate(pair = str_replace(pair, "_", " vs. "),
         average = average / overall) %>% 
  ggplot(aes(x = species, y = average)) +
  geom_boxplot(outlier.shape = NULL,
               outlier.alpha = 0,
               fill = "gray",
               alpha = 0.5) +
  geom_point(aes(fill = pair),
             pch = 21,
             color = "black",
             size = 2) +
  theme_bw() +
  coord_flip() +
  labs(x = "Spp", y = "Avg % disim")

p1

Version animada

plotly::ggplotly(p1)

Excluir localidades de Cedros en el análisis con densidades, pero incluir una perspectiva de análisis de presencia/ausencia con TODAS las localidades

Análisis de similitud

Ahora hacemos un ANOSIM incluyendo las localidades de Cedros y San Benito, pero al haber identificado diferencia sen composicion, nos quedamos solamente con presencia / ausencia.

anosim_c_cedros <- kelp %>% 
  filter(year < 2013,
         !location %in% localidades_excluidas) %>% 
  group_by(location, site, level, genusspecies) %>% 
  summarize(abundance = sum(abundance)) %>% 
  mutate(abundance = 1) %>% 
  ungroup() %>% 
  spread(genusspecies, abundance, fill = 0)

data_c_cedros_groups <- anosim_c_cedros %>% 
  select(location, site, level) %>% 
  mutate(loc_site = paste(location, site, sep = "-"))

data_c_cedros_samples <- anosim_c_cedros %>% 
  select(-c(location, site, level)) %>% 
  vegdist(method = "bray")

set.seed(43)
mds <- metaMDS(data_c_cedros_samples, trace = F)

stress <- paste("2D Stress =", formatC(mds$grstress, digits = 4, format = "f"))

cbind(data_c_cedros_groups, scores(mds)) %>% 
  ggplot(aes(x = NMDS1, y = NMDS2, sitio = site)) +
  geom_point(size = 4, aes(shape = level, color = location)) +
  coord_equal() +
  scale_color_brewer(palette = "Paired") +
  annotate(geom = "text", x = 0, y = 0.5, label = stress)

SIMPER e IVB

PCA incluyendo factores Kelp/ Temp/ Depth

FALTA INCORPORAR ESTOS DATOS, PREGUNTAR A ARTURO POR ELLOS

Pregunta 4: Que factores de hábitat explican mejor las diferencias? (Esto será mas sencillo analizar con los datos de los dos años de densidades de kelp)

Pregunta 5 Incluir ambas bases de datos 2011 y 2013

Existen cambios a traves del tiempo en esta comunidad?

Hacemos un nMDS comparandoo a nivel de año - localidad - sitio - nivel. En la primer columna se presenta la informacion total. La segunda presenta la transicion de cada localidad / sitio / nivel al conectar los puntos de de 2011 con los de 2013. A y B presentan nMDS con densidades (transformadas por raiz) y C y D presencia / ausencia. Excluimos las muestras de 2012 (Cedros) y los sitios para los que no hay muestras pareadas 2011 - 2013.

anosim_c_cedros <- kelp %>% 
  filter(location %in% c("ASA", "BMA", "ERE", "ERO", "ISME", "ITSP", "RET", "SMI", "SSI")) %>% 
  group_by(year, location, site, level, transect, genusspecies) %>% 
  summarize(abundance = sum(abundance)) %>% 
  ungroup() %>% 
  group_by(year, location, site, level, genusspecies) %>% 
  summarize(abundance = mean(abundance)) %>%
  ungroup() %>% 
  mutate(abundance = sqrt(abundance)) %>% 
  spread(genusspecies, abundance, fill = 0)

data_c_cedros_groups <- anosim_c_cedros %>% 
  select(year, location, site, level) %>% 
  mutate(loc_site = paste(location, site, sep = "-"))

data_c_cedros_samples <- anosim_c_cedros %>% 
  select(-c(year, location, site, level)) %>% 
  vegdist(method = "bray")

set.seed(43)
mds <- metaMDS(data_c_cedros_samples, trace = F)

stress <- paste("2D Stress =", formatC(mds$grstress, digits = 4, format = "f"))

p1 <- cbind(data_c_cedros_groups, scores(mds)) %>% 
  ggplot(aes(x = NMDS1, y = NMDS2)) +
  geom_point(aes(shape = level, color = location), size = 4) +
  coord_equal() +
  scale_color_brewer(palette = "Paired") +
  annotate(geom = "text", x = 0, y = 0.5, label = stress)

p2 <- cbind(data_c_cedros_groups, scores(mds)) %>% 
  ggplot(aes(x = NMDS1, y = NMDS2)) +
  geom_line(aes(group = paste(location, site, level))) +
  geom_point(aes(shape = paste(year, level), color = location, fill = location), size = 4) +
  coord_equal() +
  scale_color_brewer(palette = "Paired") +
  scale_fill_brewer(palette = "Paired") +
  annotate(geom = "text", x = 0, y = 0.5, label = stress) +
  scale_shape_manual(values = c(1, 2, 21, 24))


# Reciclamos el codigo y lo hacemos ahora por presencia / ausencia
anosim_c_cedros <- kelp %>% 
  filter(location %in% c("ASA", "BMA", "ERE", "ERO", "ISME", "ITSP", "RET", "SMI", "SSI")) %>% 
  group_by(year, location, site, level, genusspecies) %>% 
  summarize(abundance = sum(abundance)) %>% 
  ungroup() %>% 
  mutate(abundance = 1) %>% 
  spread(genusspecies, abundance, fill = 0)

data_c_cedros_groups <- anosim_c_cedros %>% 
  select(year, location, site, level) %>% 
  mutate(loc_site = paste(location, site, sep = "-"))

data_c_cedros_samples <- anosim_c_cedros %>% 
  select(-c(year, location, site, level)) %>% 
  vegdist(method = "bray")

set.seed(43)
mds <- metaMDS(data_c_cedros_samples, trace = F)

stress <- paste("2D Stress =", formatC(mds$grstress, digits = 4, format = "f"))

p3 <- cbind(data_c_cedros_groups, scores(mds)) %>% 
  ggplot(aes(x = NMDS1, y = NMDS2)) +
  geom_point(aes(shape = level, color = location), size = 4) +
  coord_equal() +
  scale_color_brewer(palette = "Paired") +
  annotate(geom = "text", x = 0, y = 0.5, label = stress)

p4 <- cbind(data_c_cedros_groups, scores(mds)) %>% 
  ggplot(aes(x = NMDS1, y = NMDS2)) +
  geom_line(aes(group = paste(location, site, level))) +
  geom_point(aes(shape = paste(year, level), color = location, fill = location), size = 4) +
  coord_equal() +
  scale_color_brewer(palette = "Paired") +
  scale_fill_brewer(palette = "Paired") +
  annotate(geom = "text", x = 0, y = 0.5, label = stress) +
  scale_shape_manual(values = c(1, 2, 21, 24))

plot_grid(p1, p2, p3, p4, labels = "AUTO", ncol = 2)

ANOSIM con años

El ANOSIM lo hacemos con densidades promediadas y transformadas por raiz cuadrada.

anosim_c_cedros <- kelp %>% 
  filter(location %in% c("ASA", "BMA", "ERE", "ERO", "ISME", "ITSP", "RET", "SMI", "SSI"),
         transect > 0) %>% 
  group_by(year, location, site, level, transect, genusspecies) %>% 
  summarize(abundance = sum(abundance)) %>% 
  ungroup() %>% 
  mutate(abundance = sqrt(abundance)) %>% 
  spread(genusspecies, abundance, fill = 0)

data_c_cedros_groups <- anosim_c_cedros %>% 
  select(year, location, site, level) %>% 
  mutate(loc_site = paste(location, site, sep = "-"),
         year_loc_site = paste(year, location, site, sep = "-"))

data_c_cedros_samples <- anosim_c_cedros %>% 
  select(-c(year, location, site, level)) %>% 
  vegdist(method = "bray")

perm <- how(nperm = 999)

adonis(formula = data_c_cedros_samples ~ year + location, data = data_c_cedros_groups, permutations = perm)

## 
## Call:
## adonis(formula = data_c_cedros_samples ~ year + location, data = data_c_cedros_groups,      permutations = perm) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##            Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## year        1     1.035 1.03547  5.3095 0.02273  0.002 ** 
## location    8     6.680 0.83497  4.2814 0.14665  0.001 ***
## Residuals 194    37.834 0.19502         0.83062           
## Total     203    45.549                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Mapa

Le falta mucho diseño al mapa, pero es la idea… hacer zoom a las areas de interes

proj <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
data(World)

World <- st_as_sf(World) %>% 
  st_transform(crs = proj) %>% 
  filter(iso_a3 == "USA")

coastline <- readRDS(here("raw_data", "spatial", "coastline_mx.rds")) %>% 
  st_as_sf()

points <- kelp %>% 
  group_by(year, location, site, latitude, longitude) %>% 
  count()

major <- ggplot() +
  geom_sf(data = World, fill = "gray") +
  geom_sf(data = coastline, fill = "gray") +
  geom_point(data = points, aes(x = longitude, y = latitude), size = 2, fill = "steelblue", shape = 21)

general <- major +
  xlim(c(-120, -110)) +
  ylim(c(27, 35))

north <- major +
  xlim(c(-117.5, -116)) +
  ylim(c(31.5, 32.5))

center <- major +
  xlim(c(-116.5, -115)) +
  ylim(c(29.5, 31.5))

south <- major +
  xlim(c(-116, -114.5)) +
  ylim(c(28, 28.5))

left <- plot_grid(general, NULL, labels = c("A", NA), ncol = 1)
right <- plot_grid(north, center, south, labels = c("B", "C", "D"), ncol = 1)

plot_grid(left, right, labels = c("A", NA), ncol = 2)

Session info

sessionInfo()

## R version 3.5.0 (2018-04-23)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 16299)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=Spanish_Mexico.1252  LC_CTYPE=Spanish_Mexico.1252   
## [3] LC_MONETARY=Spanish_Mexico.1252 LC_NUMERIC=C                   
## [5] LC_TIME=Spanish_Mexico.1252    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] bindrcpp_0.2.2     forcats_0.3.0      stringr_1.3.0     
##  [4] dplyr_0.7.4        purrr_0.2.4        readr_1.1.1       
##  [7] tidyr_0.8.0        tibble_1.4.2       tidyverse_1.2.1   
## [10] tmap_1.11-2        sf_0.6-1           cowplot_0.9.2     
## [13] ggplot2_2.2.1.9000 magrittr_1.5       here_0.1          
## [16] vegan_2.5-1        lattice_0.20-35    permute_0.9-4     
## 
## loaded via a namespace (and not attached):
##   [1] colorspace_1.3-2    deldir_0.1-15       colorRamps_2.3     
##   [4] class_7.3-14        gdalUtils_2.0.1.14  leaflet_2.0.0      
##   [7] rgdal_1.2-18        rprojroot_1.3-2     satellite_1.0.1    
##  [10] base64enc_0.1-3     dichromat_2.0-0     rstudioapi_0.7     
##  [13] lubridate_1.7.4     xml2_1.2.0          codetools_0.2-15   
##  [16] splines_3.5.0       R.methodsS3_1.7.1   mnormt_1.5-5       
##  [19] knitr_1.20          geojsonlint_0.2.0   jsonlite_1.5       
##  [22] tmaptools_1.2-4     broom_0.4.4         cluster_2.0.7-1    
##  [25] png_0.1-7           R.oo_1.22.0         rgeos_0.3-26       
##  [28] shiny_1.0.5         httr_1.3.1          compiler_3.5.0     
##  [31] backports_1.1.2     mapview_2.3.0       assertthat_0.2.0   
##  [34] Matrix_1.2-14       lazyeval_0.2.1      cli_1.0.0          
##  [37] later_0.7.1         htmltools_0.3.6     tools_3.5.0        
##  [40] coda_0.19-1         gtable_0.2.0        glue_1.2.0         
##  [43] reshape2_1.4.3      gmodels_2.16.2      V8_1.5             
##  [46] Rcpp_0.12.16        cellranger_1.1.0    raster_2.6-7       
##  [49] spdep_0.7-7         gdata_2.18.0        nlme_3.1-137       
##  [52] udunits2_0.13       iterators_1.0.9     crosstalk_1.0.0    
##  [55] psych_1.8.3.3       rvest_0.3.2         mime_0.5           
##  [58] gtools_3.5.0        XML_3.98-1.11       LearnBayes_2.15.1  
##  [61] MASS_7.3-49         scales_0.5.0.9000   hms_0.4.2          
##  [64] promises_1.0.1      parallel_3.5.0      expm_0.999-2       
##  [67] RColorBrewer_1.1-2  yaml_2.1.18         curl_3.2           
##  [70] geosphere_1.5-7     stringi_1.1.7       jsonvalidate_1.0.0 
##  [73] highr_0.6           foreach_1.4.4       e1071_1.6-8        
##  [76] boot_1.3-20         spData_0.2.8.3      rlang_0.2.0.9001   
##  [79] pkgconfig_2.0.1     bitops_1.0-6        evaluate_0.10.1    
##  [82] bindr_0.1.1         labeling_0.3        htmlwidgets_1.2    
##  [85] tidyselect_0.2.4    osmar_1.1-7         plyr_1.8.4         
##  [88] R6_2.2.2            DBI_0.8             haven_1.1.1        
##  [91] pillar_1.2.1        foreign_0.8-70      withr_2.1.2        
##  [94] mgcv_1.8-23         units_0.5-1         RCurl_1.95-4.10    
##  [97] sp_1.2-7            crayon_1.3.4        modelr_0.1.1       
## [100] rmapshaper_0.4.0    KernSmooth_2.23-15  plotly_4.7.1       
## [103] rmarkdown_1.9       readxl_1.1.0        grid_3.5.0         
## [106] data.table_1.10.4-3 digest_0.6.15       classInt_0.2-3     
## [109] webshot_0.5.0       xtable_1.8-2        httpuv_1.4.1       
## [112] R.utils_2.6.0       stats4_3.5.0        munsell_0.4.3      
## [115] viridisLite_0.3.0